RICE UNIVERSITY 


Robust Constrained Optimization Approach to Control Design 
for International Space Station Centrifuge Rotor Auto 
Balancing Control System 

By 

Barry Dirk Postma 


A THESIS SUBMITTED 
IN PARTIAL FULFILLMENT OF THE 
REQUIREMENTS FOR THE DEGREE 


Master of Science 

APPROVED, THESIS COMMITTEE: 


Pol D. Spanos, Professor, L.B. Ryon 
Endowed Chair in Engineering 
Mechanical Engineering & Material Science 


Satish Nagarajaiah, Associate Professor 
Mechanical Engineering & Material Science 


Marcia K. O’Malley, Assistant Professor 
Mechanical Engineering & Material Science 


Jiann-Woei Jang, Technical Supervisor 
The Charles Stark Draper Laboratory 


HOUSTON, TEXAS 
APRIL 2005 



Note: The views expressed in this thesis are those of the author and do not reflect the 
official policy or position of the United States Air Force, Department of Defense, or the 


U.S. Government 



ABSTRACT 


Robust Constrained Optimization Approach to Control Design for International Space 
Station Centrifuge Rotor Auto Balancing Control System 

by 

Barry Dirk Postma 

This thesis discusses application of a robust constrained optimization approach to 
control design to develop an Auto Balancing Controller (ABC) for a centrifuge rotor to 
be implemented on the International Space Station. The design goal is to minimize a 
performance objective of the system, while guaranteeing stability and proper 
performance for a range of uncertain plants. The performance objective is to minimize 
the translational response of the centrifuge rotor due to a fixed worst-case rotor 
imbalance. The robustness constraints are posed with respect to parametric uncertainty in 
the plant. The proposed approach to control design allows for both of these objectives to 
be handled within the framework of constrained optimization. The resulting controller 
achieves acceptable performance and robustness characteristics. 
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1 Introduction 

The Japanese Aerospace Exploration Agency (JAXA) plans to participate in the 
development of the International Space Station (ISS) by providing a centrifuge to be used 
for scientific experiments. The centrifuge will be housed in the Centrifuge 
Accommodation Module to be assembled with the ISS in the configuration shown in 
Figure 1-1 [1]. 



Figure 1-1. Location of the Centrifuge Accommodation Module 


The centrifuge would serve as a home for a number of rodents (or other test 
specimen). The rodents are allowed to move within habitats that are housed in the 
centrifuge rotor as shown in Figure 1-2. 
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Side View of Centrifuge Rotor Top View of Centrifuge Rotor 

Figure 1-2. Side/Top View of Rodents in Habitats 


Rodent movement within the habitats introduces an imbalance in the rotor that 
causes rotor vibration while the rotor is spinning. This imbalance must be cancelled to 
avoid the transfer of disturbances from the centrifuge rotor to the ISS. An active control 
system must be introduced to cancel imbalances in the rotor. 

Methodologies for limiting the vibration of spinning rotors have been presented in 
the literature. These methods can be divided into two techniques: direct active vibration 
control and active balancing techniques using mass redistribution. Direct active vibration 
control uses magnetic bearings that generate an external lateral force on the rotor, while 
active balancing uses a mass redistribution device to eliminate the rotor imbalance. 
Meirovitch presented an overview of active vibration control in [2], These methods can 
be applied to the direct active vibration control for a spinning rotor. Knopse presented 
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extensive work in [3]-[6] using adaptive open-loop control to suppress the rotor vibration 
using magnetic bearings as actuators. Herzog presented a method in [7] that also 
employed magnetic bearings and used notch filters in the control system. Fan introduced 
a vibration control method in [8] for an asymmetrical rotor using magnetic bearings. 

The Auto Balancing Controller (ABC) for the centrifuge rotor is classified as a 
mass redistribution type of controller. Using mass redistribution techniques, Van De 
Vegte first proposed a modal balancing method for the balancing of spinning rotors 
during operation [9]. Gosiewski proposed a control method based on modal analysis of 
the rotor vibrations using a computer as the controller in [10] and [11]. Dyer and Ni 
introduced a balancing method based on influence coefficients to achieve on-line 
estimation and active control [12]. Many of these methods were developed for the case 
of a constant rotational rotor speed. More recently, Zhou and Shi developed several 
balancing methods for speed-varying rotors presented in [13]-[15]. A thorough summary 
of work in the area of active balancing and vibration control of rotor systems is presented 
in [16]. In those methods, it is assumed that a rotor speed is either known or estimated 
exactly; the issue of control system robustness is never addressed. It is important to 
account for the linear-time varying dynamics of the rotor when using these approaches. 
In this thesis, a constrained optimization approach is proposed to design an ABC that can 
robustly cancel rotor imbalances by redistributing mass in the presence of system 
uncertainty. With the proposed methodology, both time domain and frequency domain 
stability and performance requirements can be easily included in the design process. 
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The remainder of this thesis is divided into five chapters. 

Chapter 2 provides relevant problem background. A description of the centrifuge 
rotor system and the ABC system is presented. 

Chapter 3 presents technical background for the solution of the problem. System 
robustness for multivariable systems is discussed, and the solution method is introduced. 

Chapter 4 explains the method used to solve the control problem. This section 
describes the implementation of the robust constrained optimization approach to control 
design and how it was applied to design an ABC system. 

Chapter 5 presents the results of the control design. The system is analyzed in 
both the time domain and the frequency domain. Also, the results of Monte Carlo 
simulations are presented. 

Chapter 6 summarizes the topics in this thesis and presents suggestions for 
possible future work related to this problem. 
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2 Problem Background 

In order to delineate the methodology used to solve the ABC design problem, it is 
important to properly describe the centrifuge rotor system. In this section the centrifuge 
rotor system and the ABC system are described in detail. Following the system 
description is a general statement of the problem to be solved in the design of an ABC. 

2.1 Centrifuge Rotor System 

The purpose of the centrifuge is to simulate varying levels of gravity so that the 
effects on the rodents can be studied. The centrifuge achieves the effect of varying levels 
of gravity by spinning at different steady-state rates. In order to reach a desired level of 
simulated gravity the rotor must go through a “spin-up” when the spin rate is slowly 
increased to the desired level. 

The spinning centrifuge imparts a centrifugal force on the rodents that serves to 
simulate gravity. The centrifugal force increases with increasing spin rate. An 
illustration of how the centrifugal force acts on the rodents is shown in Figure 2-1. 
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Figure 2-1. Illustration of Centrifugal Force Acting on Rodents 


The centrifugal force acting on an individual rodent is specified by the equation 

Force =M wdeM r or , (1) 

where the symbol Force refers to the centrifugal force on the rodent (N), M r0 dent is the 
mass of rodent (kg), r is the distance from spin axis to rodent (m), and co represents the 


spin rate of centrifuge rotor (rad/sec). 
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The purpose of the ABC system is to balance the rotor when the rodents have 
imposed an imbalance. An illustration of the rotor and balancer masses is given in Figure 
2 - 2 . 



Figure 2-2. Rotor and Balancer Masses 


If all of the rodents move to the same corner of their habitat, the center of mass of 
the rotor shifts in that direction. The ABC utilizes sliding balancer masses as pictured in 
Figure 2-2 within the rotor to cancel such imbalances. The static balancer masses are 
utilized to correct imbalances that cause the rotor to translate in the (£,T|) plane. The 
dynamic balancer masses are utilized to correct imbalances that cause the rotor to tilt 
about the (£,r|) axes. An illustration of the static balancer masses canceling an imbalance 
is given in Figure 2-3. The dynamic balancers work in a similar manner to cancel 


imbalances. 



1 . Balanced Rotor 



direction of the imbalance 



2. Imbalance due to 
rodent position 



4. Balancer masses move 
to correct imbalance 


Figure 2-3. Illustration of Balancer Masses Canceling an Imbalance 

The naming convention for the balancer masses and positions is described in 
Figure 2-4. There are two types of balancer masses: static balancer masses (rectangles in 
Figure 2-4) and dynamic balancer masses (circles in Figure 2-4). The set of static 
balancer masses moving in the same direction along the ^-axis is denoted by the subscript 
bl and the set of static balancer masses moving in the same direction along the rpaxis is 
denoted by the subscript b2. The symbols £bi and r|b 2 are used to describe the position of 
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the set of balancer masses along the respective axes. The dynamic balancer masses also 
work in pairs. Each mass is paired with the dynamic balancer mass diagonal to it. The 
masses in a pair of dynamic balancer masses move in opposite directions parallel to the 
^-axis (which is pointing out of the page in Figure 2-4). For instance, if the dynamic 
balancer mass in the upper right quadrant of Figure 2-4 moves up, the balancer mass in 
the lower left quadrant moves the same distance in a downward direction. These mass 
positions are denoted by £ b 3 and ^,4 corresponding to balancer mass names in Figure 2-4. 



There exists a maximum imbalance that the rodents can impose on the system. 
The maximum imbalance occurs when all of the rodents are positioned in the same corner 
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of their habitats. This moves the center of mass of the rotor the furthest distance from the 
spin axis of the rotor. Such an imbalance is illustrated in Figure 2-5. 



In this case all of the rodents are positioned in the upper-right comer of their 
respective habitats. Equally large imbalances occur if all of the rodents move to any 
other corner of their habitats at the same time. This maximum imbalance is used to 
evaluate controller performance in Chapter 4 and Chapter 5. 

Note that there are two separate coordinate frames used with respect to the 
centrifuge rotor. They are shown in Figure 2-6. The (x,y) coordinate frame is fixed with 
respect to the centrifuge base, and the (^/n) coordinate frame rotates with the centrifuge 


11 


rotor. Both coordinate frames are used since the rotor displacement sensors work in the 
fixed coordinate frame and the input to the ABC is a measurement in the rotating 
coordinate frame. 



Figure 2-6. Rotor Coordinate Frames 


A coordinate transformation as shown in Equation ( 2 ) is used to change 
coordinate frames. Specifically, 

The ABC system is used to limit the translational response of the rotor due to 
imbalances imposed by the position of the rodents within the habitats. The input to the 
controller is an absolute measurement of the rotor displacement (translation and tilt) as 
shown in Figure 2-7. The centrifuge rotor sensors measure relative displacement 
between the rotor and the International Space Station. With this relative measurement. 


£\ cos(ryr) sin(ryr)] x\ 

r}\ _-sin(tyr) cos(<2*)^|_ 


( 2 ) 
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the absolute rotor displacement can be estimated using a Kalman filtering technique [17]. 
In this thesis it is assumed that these absolute rotor measurements are already available. 



Figure 2-7. Rotor Translation and Tilt 


The equations of motion for the centrifuge rotor in the fixed frame coordinate 
system (x,y coordinate system) are derived below. First, the equations of motion in the 
translational direction are derived for the system in Figure 2-8. Then, the equations of 
motion in the tilting direction are derived for the system in Figure 2-9. 
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Figure 2-8. Model for Derivation of Equations of Motion in Translational Direction 


The equations of motion in the translational direction [17] can be written in each 


axis as 


F x = Mix 
and 

F y = My. 


(3) 

(4) 


The symbol F denotes the external force on the rotor due to springs, dampers, 
control inputs, and external disturbances. The equations of motion can thus be written 

F x = -C x (x-Lsin<f> y )- K x (x-Lsm<p y ) + u x +d x ( 5 ) 


and 


F y = —C y (y + Lsin0 x ) — K y (y + Lsin 0 X ) + u y +d y . 


( 6 ) 
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where the symbols x and y denote the translation of the rotor in the fixed coordinate 
frame, <|> x and <|> y are the tilt of the rotor in the fixed coordinate frame, u is the input force 
due to balancer mass position, d is the disturbance force due to rodent position, M is the 
mass of the rotor, C and K denote damping and stiffness coefficients between the rotor 
and the outer wall of centrifuge, and L is the distance from the rotor center of mass to the 
base of the rotor where the springs and dampers act. 

Using the small angle approximation and making the substitutions 

C„=C x =C, t (7) 

K u =K x =K y , (8) 

C„ =-C x L = -C y L, (9) 

and 

K n = -K x L = -K y L, (10) 

the equations of motion in the translational direction can be written in the form 

Mx + C n x + C n <j> y + K u x + K n (j) s =u x +d x ( 11 ) 

and 

My + C n y-C n <j> x + K n y- K n 0 x =u y +d y . ( 12) 

To derive the equations of motion in the tilting direction, the model in Figure 2-9 
is considered. The symbols Kry and C ry denote rotational stiffness and rotational 
damping, respectively, about the y-axis. The translational springs/dampers in the y-axis 
and the rotational springs/dampers about the x-axis are not shown in the figure. 
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Figure 2-9. Model for Derivation of Equations of Motion in Tilting Direction 


With the small angle approximation, the equations of motion in the tilting 
direction can be written in the form [19] 


IJ V -I z co<j> x 


= T , 


(13) 


where 0) is the spin rate of the rotor, la is the transverse moment of inertia of the rotor, I z 
is the polar moment of inertia of the rotor, and T is the external torque applied to the 


rotor. 
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Specifically, the symbol T denotes the external torque on the rotor due to springs, 
dampers, control inputs, and external disturbances. The torque can be written as 

T * = ~ C y L y - C rA - C y L 2 sin <p x - K y Ly - K rx <f> x - K y L 2 sin </> x +u^+d^ ( 14 ) 

and 

= C x Lx - Cj y - C X L 2 sin <p y + K x Lx - K ry <p x - K X L 2 sin <p y +u^+d^, (15) 

where u^, is the input torque in the indicated direction due to balancer mass position and 
d^, is the disturbance torque in the indicated direction due to rodent position. 

Using the small angle approximation the torque can be written as 

T + =-C y Ly-Cj x -C y L 2 ^-K y Ly-K rx <P x -K x L 2 <P x +u^+d^ ( 16) 

and 

T„ - C,U-Cj, -C.L% * - K,L% + «*+<**■ ( 17 ) 

Make the substitutions 

C 22 = C n + C y L = C 0 . + C X L , (18) 

K 22 =K rx +K y L 2 =K iy +K x L 2 , (19) 

C n =-C x L = -C y L, (20) 

and 

K n =-K x L = -K y L. 


( 21 ) 
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Then, the equations of motion in the tilting direction can be written in the form 


I d 0 x ~ C \2 y + C 220x + 0)1 z 0 y ~ K ny + K 220x = U <fc + d <t* 


( 22 ) 


and 


I j(j) y + C\ 2 x + C 22 0 y col 7 <j) x + K i2 x + K 22 (f) y — u ^ + d ^ . 


(23) 


Now, the full equations of motion can be written 


Mix + Cjjir + C n 0 y + K\^x + K n 0 y — u x + d x , 


(24) 


My + C u y~C n <p x + K u y- K n 0 x =u y +d y . 


(25) 


IJ X C n y + C 22 (f) x + OJl z 0 y K ]2 y + K 22 0 x — u ^ + d ^ , 


(26) 


and 


l Ay ■*" ^ 220 y ^ z0x K] 2 X+ K 22 0 y — + d ^ . 


(27) 


The equations of motion in the rotating coordinate frame (see Figure 2-6) can be 
derived from the equations of motion in the fixed frame given in Equations ( 24 ), ( 25 ), 
( 26 ), and ( 27 ). The importance of deriving the equations of motion in the rotating 
coordinate frame is emphasized in Section 2.2. 

Begin with the equations of motion in the fixed frame (Equations ( 24 ), ( 25 ), 
( 26 ), and ( 27 )) and apply the coordinate transformations 


_yj Wj 


Y| |Y1 

I =R I 


(28) 


and 



( 29 ) 
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where 


cos(ryf) -sin(&*)l 
sin(ryr) cos(ax) 


Also, introduce R in the transformations 


and 


where 


and 



kj. 


R=RI = 


sin(ryf) 

-cos(ryf) 


cos(ryr)! 

sin(ryr)[j 


I = 




( 30 ) 


(31) 


(32) 


(33) 


(34) 


Differentiation in time of Equation ( 28 ) and Equation ( 29 ) yields 


>1 

J] 


= R 


Jf\ 


+ R 


and 


Aj 


= R 


Vj 


+ R 


J1 J 

Aj. 


(35) 


(36) 


19 


Differentiation in time of Equation ( 35 ) and Equation ( 36 ) yields 


and 


'i 1 

.n 


= R 


jr\ 


+ 2 R 


'<ij 


+ R 


' 5 ! 


■&1 •• 
x =R 


&,1 


+ 2R 


+,\ L^J KJ 


I 


+ R 


4>c "I 

AJ. 


Differentiation in time of R from Equation ( 30 ) yields 


R = 


-rysin(ryr) -tycos(tyf)] 
rycos(ryr) -rysin(ryr)|j 


= -coRI 


Differentiating in time of Equation ( 39 ) yields 


R = 


- co 1 cos(tyr) or sin(ryr) 1 _ 
-ty 2 sin(ryr) -or cos(ryrJJ 


Also, R is defined in the same manner as R where 


R=RI =-{coRI)I = (oR = 


O)COS(OX) 

G)sin(ox) 


- rysin(tyr)] 
<ycos(ruf) 


(37) 


(38) 


(39) 


(40) 


(41) 
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Then, Equations ( 24 )-( 27 ) can be recast in the form 


and 


r v 


M 


R 


+ K n R 


m 

W\ 


+ 2 R 


+ K^R 


^ + R 

W 


<?' 

VjJ 


/ r 


+ Ci 


/? 


L?j IVj 


+ c 


12 


/? 


A.1 - 
rc +R 

IAJ 


lV. 


= /? 


M/- + J 


u n + VJ 


r 

/? 


0,1 

.f +2R 

t\ L^J 


0fl . fi 


+ /? 


|Vj 


+ c 


12 




WJ 


/7j 


+ C 


22 


/? 



1 . 
+ R 

1 

+ col. 

f 

R 

Kt> 

• 

1 

1 

+ R 

i 

+ K\ 2 

f 

-R 

rcn 


J 

yn\j 

z 

\ 

IA 

J 

yn J J 


K 

W\) 


*+* K 22^ 


W\ 


= R 


U dC + d,if ,1 


1 <K ' 


With a similarity transformation [20] 


and the identities 


R~ l 


cos (cot) sin(tyr)! 

-sin(ryr) cos(&»)^J 


R~'R = 


01 


R'R = R-'(-coRl) = -6jR-'RI = -col = 


0 -co 1 

.» o'j. 


R'R =R-'(-co 2 R)=co 2 R~'R = 


-co 2 0 1 
0 -cor\ 


R-'R =R'R1 =1 = 


-° ,i 
-1 o'j 


J (42) 

(43) 

(44) 

(45) 

(46) 

(47) 


(48) 


21 


and 


R-'R =R-'{o)R)=coR~'R = 


co 0] 

o Jj 


Equations ( 42 ) and ( 43 ) can be written in the form 


( 

M 


0 -2co\ 

f 


'-co 2 

°1 

Yll 

f 

+ c, , 

W + 

'0 



\ 

Uj 

2<y 0 'j 

Wi 

0 

-A 

A 

V 

Uj 

co 

o'j 

WjJ 


+ C 


12 


T° V 

i 

• 2^- 

1 

H- 

<y 0] 
| 

S’ 

L 

-h ATj j 

^+^ 12 

r° ii 

WJ 

J.- 1 <>j 


J 

0 aJj 

l*i JJ 


Wj 

-1 o'j 

W\ L 


Ur + d r | 


and 


0,1 


+ 


0 - 2co] 


0,1 


<t> n \ L 2 *> 0 JKj L 0 -*> 2 jKj 


-co 


0 J0J 


+ c 
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Collecting terms, the final coupled equations in the rotating frame become 
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If the coupling terms are neglected, the equations become 
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where ^ and r| represent translation of the rotor in the rotating coordinate frame, 
(|); and <J)r, represent tilt of the rotor in the rotating coordinate frame, u is the input 
force/torque due to balancer mass positions (Cbi 41 b 2 ,£b 3 ,£b 4 ), and d is the disturbance 
force/torque due to rodent position. 

Note that the dynamic behavior depends on the spin rate of the rotor, to. Since a) 
can vary in time, the equations of motion represent a linear time-varying system. This 
requires that any ABC must stabilize the system for the range 0 < to < GW. 
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It can be verified that the system with uncoupled dynamics (Equations ( 54 ) and 
( 55 )) closely approximates the system with coupled dynamics (Equations ( 52 ) and 
( 53 )) by comparing the frequency response of each system. A comparison of the 
frequency response of each system in one of the translational directions (Q and one of the 
tilting directions (<|)^) is presented in Figure 2-10 and Figure 2-11. 



Figure 2-10. Bode Comparison - Coupled vs. Uncoupled Translational Direction 
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It is clear that the frequency response of the uncoupled system closely matches 
the frequency response of the coupled system. This verifies that the uncoupled system 
can be used to approximate the coupled system. This approximation is utilized in the 
controller design presented in Chapter 4. 

2.2 Auto Balancing Controller 

As previously mentioned, the objective of the ABC system is to limit the 
translation of the rotor due to imbalances imposed by the position of the rodents within 
the habitats. This includes limiting both the peak translation of the rotor and the steady- 
state translation of the rotor due to a fixed imbalance. Additionally, the control system 
should be robust to uncertainty in the rotor spin rate, co, over the range 0 < © < aw, and 
be robust to other forms of uncertainty in the plant. Finally, the controller should avoid 
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commanding the balancer masses to change position quickly, causing saturation in the 
balancer mass speed. By avoiding actuator saturation, the linearity of the control system 
is preserved. 

Based on the rotor displacement, the ABC commands the position of the balancer 
masses in order to cancel a rodent imbalance and limit the translation of the rotor. A 
simplified block diagram of the linear system is shown in Figure 2-12. 



Figure 2-12. Control System Block Diagram 


As discussed in Section 2.1, the ABC system operates in the rotational coordinate 
system. Since the absolute displacement is estimated in the fixed coordinate system [17], 
a coordinate transformation as shown in Equation ( 2 ) is required convert the sensor 
displacement to the rotating coordinate frame before it is used by the ABC. For the 
design of the ABC system, the plant dynamics can be formulated in the rotating 
coordinate frame yielding Equation ( 52 ) and Equation ( 53 ). While the ABC can be 
designed using the dynamics in the rotating coordinate frame, the implementation of the 
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controller is performed on the fixed frame dynamic model with the coordinate 
transformation. 

It is necessary to design the controller using the rotating coordinate frame to 
effectively cancel imbalances and limit the steady-state translation of the rotor. To 
illustrate this, it is necessary to see a plot of the rotor translation in each coordinate frame. 
Plots of the rotor translation in each coordinate system are shown in Figure 2-13 for a 
stationary rodent imbalance (as described in Figure 2-14) in the open-loop (uncontrolled) 
system. The rotor is rotating at a constant spin rate. 



Figure 2-13. Comparison of Rotor Translation in Fixed and Rotating Coordinate 
Frames 


It is worth noting that in the fixed frame coordinate system (x-direction), the 
response is oscillatory, while in the rotating frame coordinate system (^-direction) the 
response is nearly constant after the transient effects subside. This is a result of the fact 
that the rotating coordinate system is rotating with the rotor imbalance. The rotor 
translates in the direction of the imbalance which remains the same in the rotating frame. 
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but changes in the fixed frame. This effect is illustrated in Figure 2-14 where the 
imbalance is initially (at time = tl) in both the x-direction and the ^-direction, but as the 
rotor rotates through angle co(t2-tl) the imbalance is no longer in the x-direction but it 
remains in the ^-direction. 


y/n y 




Figure 2-14. Comparison of Imbalance in Fixed and Rotating Coordinate Frames 


The importance of this result is that for steady-state error to be limited by the 
control system, the integral of the rotor displacement must increase as long as the error 
remains present. Since the response in the fixed frame coordinate system oscillates about 
0, the integral is not increasing and the steady-state error cannot be eliminated. The 
integral of each signal is shown in Figure 2-15. Clearly, the rotating coordinate system 
must be used to limit steady-state error since the integral of the displacement is not 
increasing in the fixed frame coordinate system. 
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Figure 2-15. Comparison of the Integral of Rotor Translation in Fixed and Rotating 
Coordinate Frames 


Figure 2-16 is a block diagram of the controller. The symbols Kp, Kd, and Ki 
represent constant gain matrices of size 4x4. 



Figure 2-16. Block Diagram of Controller 

A modified integrator is used, as opposed to a true integrator (1/s), in order to 
avoid adding poles on the jco-axis of the complex plane. By adding the poles slightly in 
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the left-half-plane, the robust stability of the control system is improved. An example of 
this improvement is presented below. 

To illustrate the improvement in robust stability, a single-input, single-output 
(SISO) system is considered with pure integral control (K/s) and modified integral 
control (K/(s+p)). Block diagrams of the system in each configuration are shown in 
Figure 2-17 and Figure 2-18. 




Figure 2-18. System with Modified Integral Control 
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The symbol P denotes the plant and is specified by the equation 



m m 


The system parameters are assumed to be K=l, p=l, m=l, c=2, and k=l. A plot 
of the poles of the transfer function from r to y for each system is given in Figure 2-19. 


Poles for Pure Integrator and Modifed Integrator 



Real Axis 

Figure 2-19. Poles for Systems with a Pure Integrator and a Modified Integrator 


It is observed from Figure 2-19 that the system with the modified integrator has 
closed-loop poles farther from the j co-axis than the system with the pure integrator. In 
general, this corresponds to a system with better robust stability since it will take a larger 
perturbation for the closed-loop poles to move to the unstable region on the right side of 
the jco-axis. This is demonstrated by allowing the damping coefficient to decrease until 
instability is reached. For the system with the pure integrator the pole plot of the closed- 
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loop system is presented in Figure 2-20. The plot shows the poles for decreasing values 
of the damping coefficient, c. 



Figure 2-20. Poles for Pure Integrator with Perturbation in c 

From Figure 2-20 it is seen that the system poles cross into the unstable region for 
c = 1.0. Therefore, the pure integral control is able to tolerate a 50% variation in the 
damping coefficient. 

For the system with the modified integrator the pole plot is presented in Figure 
2-21 for decreasing values of the damping coefficient, c. 
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Figure 2-21. Poles for Modified Integrator with Perturbation in c 

From Figure 2-21 it is observed that the system poles cross into the unstable 
region for c = 0.4. Therefore, the modified integral control is able to tolerate an 80% 
variation in the damping coefficient. This is significantly larger than the system with 
pure integral control, indicating that the system with modified integral control is more 
robust to variations in the damping coefficient (as well as other plant parameters). This 
increase in robust stability comes at the cost of increased steady-state error as discussed 
in Section 3.5. 

Block diagrams of the centrifuge rotor control system are shown in Figure 2-22 
and Figure 2-23. Figure 2-22 is the system used to design the controller. The dynamic 
equations of motion are in the rotating coordinate frame and time delays are neglected. 
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Since the time delays are small, they are neglected in the linear stability analysis, but they 
are accounted for in the time domain simulation portion of the design in Chapter 4. 



Controller 


Figure 2-22. Block Diagram of System used for Design 


Figure 2-23 shows the system used to implement and test the controller. The 
system is simulated using a solver based on the Runge-Kutta (4,5) numerical integration 
scheme [21]. The dynamic equations of motion are in the fixed coordinate frame and 
time delays are present. A coordinate transform is used to transit from the fixed 
coordinate frame to the rotating coordinate frame. The symbols (Cbi.T|b 2 ,£b 3 »£b 4 ) denote 
the commanded positions of the balancer masses. 
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Figure 2-23. Block Diagram of System used for Implementation 


2.3 General Problem Statement 

It is assumed that the following constraints apply to the ABC: 

I. limit the peak translation ( 2 +tj 2 ) of the rotor < yi for all rodent 
disturbances, 

II. ensure stability over the range 0 < (0 < ow, 

III. be robust in the presence of spin rate(w), rotor mass (M), and rotor 
inertia(I<j,Iz) uncertainty, 

IV. limit commanded balancer mass speed < 72, and 

V. limit steady-state rotor translation < 73 for a fixed disturbance (rodents are 
stationary), 

where the symbols 71,72, and 73 represent specific problem constraints to be defined in 
Section 4 . 1 . Note that the tilt of the rotor is not directly constrained since it contributes to 


the translation of the rotor (which is constrained). 
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3 Technical Background 

Before introducing the design process of Chapter 4, some background in the 
technical tools used to analyze the control system is necessary. In this section, 
background in robust control for multi-input, multi-output (MIMO) systems is presented. 
Also, steady-state error in control systems is examined. Finally, a robust constrained 
optimization approach to control design is presented. 

3.1 Multi-Input, Multi-Output (MIMO) Control 

It is well known from classical control, single-input, single-output (SISO) 
systems, that the gain of a system can be defined in terms of the frequency response or 
Bode plot [22]. Defining the gain of a system for the MIMO case is less clear since a 
Bode plot only represents the frequency response from a single input to a single output. 

To describe system gain for the MIMO case, singular values are introduced [23]. 
Consider the transfer function matrix G with m inputs and n outputs. If G is evaluated at 
a fixed frequency, to, G(jco) is a constant n x m complex matrix. Any matrix G can be 
decomposed into its singular value decomposition (SVD) [23] 

H 

G = imv t ( 57 ) 

where the symbol E is an n x m matrix with non-negative singular values, G;, in 
descending order along the main diagonal, U is an n x n unitary matrix of output singular 
vectors, and V is an m x m unitary matrix of input singular vectors. 

The column vectors of U, uj, represent the output directions of the plant, while the 
column vectors of V, Vj ; represent the input directions of the plant. That is, if an input in 
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the direction of Vi is considered, the output will be in the direction u, and the signals are 
related by the corresponding singular value. Specifically, 

GVj = <7 i u i . ( 58 ) 

It can be shown that the largest gain in any input direction is equal to the largest 
singular value [24]. The largest singular value of a matrix is denoted by the symbol 

cf(G) = <7,(G). (59) 

To obtain the MIMO frequency response for G(s), G is evaluated at s=jG> for the 
interval © for which the frequency response is desired. The maximum singular values are 
then plotted as functions of © to provide the maximum system gain at any frequency, ©. 
The peak of this frequency response plot is known as the H-infinity norm of the 
system[24], defined by the equation 

||G(s)L = max5 f (G(yfl>» (60) 
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To illustrate, the system [25] in Figure 3-1 is examined. 



Figure 3-1. Control System Example 


In this example, H and K are specified by the equations 

9 -10 1 

5 + 1 5 + 1 

-8 9 

5 + 2 5 + 2 J 

and 





The transfer function from r to y can be written in the form 

y = Tr, 


where 


(61) 


(62) 


(63) 


T = (I + HK)~ l H . 


(64) 
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This matrix is evaluated at s=jco=0 to find the singular values of the system at 
co=0 rad/sec. This yields 


T(ja> = 0) 


0.6452 -0.6774] 

-0.2903 0.3548 [j 


(65) 


Using SVD, this matrix can be written in the form 


' 0.6452 -0.6774] 

'- 0.8983 

0.4395] 

"1.0413 

,1 

"- 0.6790 

0.7341] 

-0.2903 0.3548 

_ 0.4395 

0.8983^ 


0.0310] 

_ 0.7341 

0.6790] 


It is seen that the singular values at co=0 are Gi=1.0413 and G2=0.0310. Using 
SVD for each frequency in a set of desired frequencies and plotting the maximum 
singular value produces the frequency response plot in Figure 3-2. 


Singular Values of T 



Figure 3-2. Maximum Singular Values 
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It is verified that the maximum singular value at co=0 is Oi=1.0413. Singular 
values have two distinct advantages over eigenvalues for MIMO systems; they provide 
information about the gains of the plant and they can be obtained for non-square plants. 

3.2 Uncertain Systems 

One of the problem constraints presented in Section 2.3 was that the system with 
the ABC must be robust to uncertainty in plant parameters. Uncertainty in a control 
system is the discrepancy between the actual system and the model of the system. 
Uncertainty can arise from many sources. Some sources of uncertainty are [25]: 

I. approximate or erroneously valued parameters in the plant, 

II. ignored complexities (nonlinearities, time delays, etc.), 

III. imperfections in sensors/actuators, 

IV. unmodeled dynamics, and 

V. model order reduction of the plant. 

In general, uncertainty appears in two basic forms: parametric uncertainty and 
unmodeled dynamics [28]. Parametric uncertainty will be the form of uncertainty 
primarily examined in the design process for this thesis. 

Uncertainty can be easily modeled in a feedback manner known as a linear 
fractional transformation (LFT) [24] [28] [29]. Using LFTs it is possible to separate what 
is known (model of the system (M)) from what is unknown (the bounded uncertainty (A)) 
in a feedback interconnection [28] as shown in Figure 3-3. 
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Figure 3-3. Linear Fractional Transformation 

In the case of unmodeled dynamics, A is any fully populated transfer function 
matrix where 

ML SI. (67) 


This is known as unstructured uncertainty. 

In the case of parametric uncertainty, A generally has a diagonal structure of the 

form 



( 68 ) 


Equation ( 67 ) still applies to the case of parametric uncertainty, but the structure of the 
uncertainty is constrained to the form in Equation ( 68 ). This is known as structured 
uncertainty. Generally, weighting functions are added to the plant so that Equation ( 67 ) 


is satisfied. 
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To illustrate how system uncertainty can be formulated in an LFT framework, the 
system in Figure 3-4 is considered. This is the same system considered in Section 3.1 
with an uncertainty added in the input to the plant (or the actuator signal). 



Figure 3-4. Uncertain System Example 


The symbols H and K are as defined in Section 3.1 and A has a diagonal structure. 
Specifically, 



(69) 


and 5i and 62 are allowed to vary on the interval [- 0 . 8 , 0 . 8 ], allowing for 80% uncertainty 
in each actuator signal. The uncertainties can be “pulled out” and the system and can be 
re-arranged as in Figure 3-5. 
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Figure 3-5. Pulling out the Uncertainty 


The symbol M represents the transfer function from ua to yA. Using block 
diagram algebra, an expression for M is derived. Specifically, 


> 

II 

1 

(70) 

y A =~KH(u A + y A ) , 

(71) 

y A =-KHu A -KHy A , 

(72) 

(/ + KH)y A = —KHu a , 

(73) 

y A =-(I + KH)- i KHu A , 

(74) 


M = -(/ + KH)-'KH . 


(75) 
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Now, the system can be written as an LFT as in Figure 3-3. However, A must still be 
normalized so that the norm bound on A is 1. Weighting functions, Wi and W 2 , as shown 
in Figure 3-6, are used to accomplish this scaling. 



Figure 3-6. Normalizing the Uncertainty 


In this case, the weighting functions 


and 


0.8 0 1 

f 

0 0.8J 


W 2 = 


1 I 1 

L° i'j 


(76) 


(77) 


normalize the uncertainty while allowing the actuator gain 80% uncertainty. The 
weighting functions are appended to M to recast the system in the form of Figure 3-3 


with A normalized to 1. 
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3.3 Small Gain Theorem 

The small gain theorem ties together the concepts of singular values and system 
uncertainty to form the basis of the robust stability problem for MIMO systems. Robust 
stability is achieved when all plants in a given uncertainty set are stable for all possible 
perturbations [24]. Consider the system shown in Figure 3-7 where M(s) and A(s) are 
strictly stable transfer function matrices. 



The small gain theorem [29] states that for y>0 the interconnected system in 
Figure 3-7 is well-posed and internally stable for all A(s) with 

||A|L ^ 1 1 Y tf an d only if ||M|| ot <y. ( 78 ) 

This condition can be re-written as 


c t(M (jco))cr(A(jo))) < 1 V<w. 


(79) 
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Commonly, the problem is posed such that A is norm-bounded by 1. In this case, 
the robust stability (RS) condition based on the small gain theorem can be written as 

RS a{M{jco)) < 1 Vo. ( 80 ) 

Alternatively, the condition can be written as 

RS<=>||mL<1. (81) 

3.4 The Structured Singular Value 

The small gain theorem provides the foundation for answering the problem of 
robust stability. However, the small gain theorem is conservative for the case of 
structured uncertainty, or parametric uncertainty, since it assumes an unstructured 
perturbation. The structured singular value is an expansion of the small gain theorem 
used to provide a less conservative answer to the question of robust stability (or robust 
performance) subject to structured uncertainty. The structured singular value was first 
introduced by Doyle in 1982 [26] [27]. 

3.4.1 Robust Stability with the Structured Singular Value 

The small gain theorem gives a condition for robust stability in Equation ( 80 ). 
While this condition provides a method for evaluating robust stability, it is known to be 
conservative for the case of structured uncertainty [25]. The structured singular value (|i) 
is introduced to provide a less conservative answer for the case of structured uncertainty. 
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The structured singular value, jll, is defined as 

^ ^ min{a(A) s.t. det (/ - MA) = 0} ’ ^ ^ 

where the denominator is the maximum singular value of the smallest A that is able to 
destabilize the system. While it is quite difficult to compute p directly, advanced 
algorithms exist [28] to compute upper and lower bounds on p. This thesis will be 
primarily concerned with the upper bound since it provides a slightly conservative 
answer, while the lower bound often provides an answer that is not conservative enough. 
From the small gain theorem, an initial upper bound on p can be formulated. 
Specifically, 

ju A (M)<a(M(ja»). (83) 

As previously stated, the bound in Equation ( 83 ) is known to be conservative for 
the case of structured uncertainty. To reduce conservatism, scaling matrices, D and D 1 , 
are introduced as shown in Figure 3-8. Figure 3-7 is replaced with Figure 3-8 for 


evaluating robust stability with p. 
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A’ 


M’ 

Figure 3-8. Robust Stability with Scaling Matrices 

Since A has a diagonal structure for the case of structured uncertainty and D is a block 
diagonal scaling matrix, 

DAD' 1 = A . ( 84 ) 

The block diagram can now be drawn as in Figure 3-9. 


M’ 

Figure 3-9. Simplified Robnst Stability with Scaling 
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The small gain theorem now provides the upper bound on p 

ju(M) < (j{DMD ~ x ) (85) 

where D is a free parameter that can be used to get a less conservative upper bound on p. 
It follows that an optimal D can be found to find the least conservative upper bound on p, 
which provides a new, less conservative condition for robust stability. That is, 

RS o min DMD~ X <1 ( 86 ) 

D 

To illustrate the advantage of using the structured singular value to determine 
robust stability, the example from Section 3.2 is again considered. The system with 
weighting functions as specified in Section 3.2 is presented in Figure 3-10. The system, 
M, is taken to be the closed-loop transfer function, T, as defined in Section 3.1. 



Figure 3-10. System with Weighting Functions 
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For evaluating robust stability, there are two conditions. Both provide an upper 
bound on |i. Specifically, 

RS <=> cr(M(jot)))< 1 Vty Small Gain Theorem ( 87 ) 


and 


RS <=> min 

D 


DMD~ X 


<1 


Structured Singular Value 


( 88 ) 


The upper bounds on |i for each case are plotted in Figure 3-11. In this case the 
small gain theorem does not guarantee robust stability for the system since the peak value 
is greater than 1. However, the structured singular value does guarantee robust stability 
since the peak is below 1. Clearly, the small gain theorem is conservative in this case. 
The structured singular value allows for a less conservative answer to the question of 
robust stability through the use of the scaling matrices, D and D 1 . 


Small Gain Theorem vs. p. 
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3.4.2 Robust Performance with the Structured Singular Value 

Robust performance is achieved when a given performance objective is achieved 
for all plants in a given uncertainty set. More specifically, for the system in Figure 3-12, 
robust performance is defined as 

RP o ||FL < a . ( 89 ) 



Figure 3-12. Definition of Robust Performance 


Generally, weighting functions are added to F so that the problem can be written 

RP «=> |f|L < 1. (90) 

The problem of robust performance can be treated in the framework of the 
structured singular value by using the small gain theorem to convert the problem to a 
robust stability problem. Since the requirement for robust performance in Equation ( 90 ) 
is the same as the small gain requirement for robust stability in Equation ( 81 ) where M 
is replaced by F, the robust performance problem can be posed as a robust stability 
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problem. Figure 3-13 is a pictorial representation of a robust performance problem being 
converted to a robust stability problem through the use of a fictitious uncertainty block A P 
[24]. 



Figure 3-13. Converting a RP Problem to a RS Problem 


Now, if the robust stability requirement 

RS^> J u l (M)<l \/Q) (91) 

is satisfied, robust performance is guaranteed. To illustrate how robust performance is 
evaluated with the structured singular value, the example from Section 3.2 is again 
examined. The system in Figure 3-14 is considered. 
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The symbol n represents the system input, sensor noise, and the symbol y 
represents the variable to be regulated. It is assumed that for the system with 80% 
actuator uncertainty, the performance objective is to limit the output, y, to 10% of the 
sensor noise, n, at high frequencies (where sensor noise is most likely). This 
performance objective requires that the system is robust to the sensor noise. To 
normalize the performance objective to 1, weighting functions are added to the system as 
shown in Figure 3-15. 
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t 



Figure 3-15. System with Weighting Functions 


The symbols Wi and W 2 are as described in Section 3.2. The matrix 




(92) 


serves to impose the 10% limit on the output of the system. Further, W n serves to impose 
the robust performance requirement at high frequencies only. In this case the output is 
required to be less than 10% of the sensor noise at frequencies above 100 rad/sec 
(roughly where the plot of W„’ (shown in Figure 3-16) nears 1). Specifically, 


W = 

r n 


W. 


1 

wjj 


(93) 


and 


W '= 

n 


5 + 1.73 
5 + 173 . 


( 94 ) 




Magnitude (abs) 
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Weighting Function Wn' 



Figure 3-16. Weighting Function W„' 


This performance objective is evaluated using p. The upper bound of p is plotted 


for the 10% requirement in Figure 3-17. 
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Robust Perfomance - Limit error to 1 0% 



Figure 3-17. Robust Performance Plot for 10% Requirement 


Since the peak of the robust performance plot is greater than 1, the conditions for 
robust performance are not satisfied. However, if the performance requirement is relaxed 
to limit the output to 20% of the input signal over the same frequency range, robust 
performance is achieved. The results are shown in Figure 3-18. 



56 


Robust Performance 



Figure 3-18. Comparison of 10% and 20% Requirements for Robust Performance 


3.4.3 Real Parametric Uncertainty to Complex Uncertainty 

For a complex uncertainty (unmodeled dynamics or complex parameter 
uncertainty) p is a continuous function. However, p is not necessarily a continuous 
function when all perturbation blocks are real, as is often the case for parametric 
uncertainty [30). This discontinuity can cause problems in the convergence of the bounds 
on p [28]. 

In fact, in order to derive a more reliable bound for the design of the ABC, the 
parametric uncertainty in the spin rate was treated as an uncertain complex parameter. 
The interval for spin rate was then divided into several smaller complex intervals and 
each interval was tested for robust stability and performance [31]. An illustration of how 
the uncertainty in the spin rate was treated is shown in Figure 3-19. 
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Real Perturbation in co Complex Perturbations in co 



Figure 3-19. Real to Complex Uncertainty 


The complex uncertainty representation is somewhat conservative since it allows 
for a small complex contribution in the spin rate, which is not possible. However, the 
complex representation contains the entire interval in the real uncertainty representation 
and it allows for a more reliable computation of the bound on the structured singular 
value. 

3.4.4 Centrifuge Rotor System in |i-analysis Framework 

Figure 3-20 presents a pictorial representation of an uncertain centrifuge rotor 
system formulated into a robust stability problem using the structured singular value. 
The diagram describes the steps necessary to put the uncertain centrifuge rotor system in 
the p-analysis framework. 
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Now, the plant dynamics 
can be put into the 
following form 


&{G(jco))cr(A(jG)))<l Vry 

Structured 
Singular 
Value (|i) 4 , 

ti{G(jaj))a{A(jcoj)<\ Vco 



If A is normalized to 1 , 
the requirement for 
robust stability is: 


H(G(jco)) < 1 


P = Plant 
K = Controller 


Figure 3-20. Evaluating System Robustness with the Structured Singular Value 


3.5 Steady-State Error in Control Systems 

One of the problem constraints as presented in Section 2.3 is to limit the steady- 
state error in the centrifuge rotor system. It is common practice in control design to 
introduce an integral term into the controller to eliminate steady-state error. This integral 
term increases the system type. The system type is defined as the number of poles the 
open-loop system (the transfer function from e to y in Figure 3-21) has at s = 0 [22]. 
Adding an integral term in the controller introduces a single pole at s = 0. The steady- 
state error of a system can be expressed with respect to the system type and the input to 
the system. For example, the steady-state error due to a step input for a type 0 system 
and type 1 systems and higher is shown in Table 3-1. 
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System Type 

Steady-State Error 

Type 0 System 

e ss = constant 

Type 1 (or higher) System 

Css = 0 


Table 3-1. Steady-State Error due to Step Input 


A controller utilizing only an integral term has the structure shown in Figure 3-21. . 
P is a type 0 system. Since an integral term has been added via the controller to a type 0 
system, the system is now a type 1 system. It should have steady-state error of 0 for a 
step input. 



Figure 3-21. System with Pure Integral Control 


If K is chosen so that the system is stable, the integral controller in Figure 3-21 
will eliminate steady-state error. However, adding an integral term to the controller adds 
a pole of the open-loop system on the j co-axis. Although this pole can be moved off the 
axis for the closed-loop system, K is often constrained to be small to satisfy other 
problem constraints. In this case, the pole is not moved far off of the j co-axis and this 
close proximity to the unstable region can lead to poor system robustness. For this 
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reason, it is often advantageous to use a modified integrator, such as in Figure 3-22, to 
move the pole away from the jo-axis when p is a positive real number. The system is a 
type 0 system and steady-state error is expected for a step input. 



The increase in system robustness for modified integral control comes at the cost 
of allowing some steady-state error for a step input. This relationship is derived through 
use of the final value theorem [22] which provides the relationship 

lim/(r) = lim^F(i) (95) 

f— S— »0 


The final value theorem relates the final value, or steady-state value in the time 
domain, to the Laplace transform of a system evaluated at s=0. As an example, consider 
the system presented in Figure 3-22 where the plant is defined by the equation 


P = 


2 C 
s + — s + — 

m 


k 

m 


(96) 


Since the steady-state error is of interest, the transfer function from r to e will be 
evaluated. The closed-loop transfer function from r to e can be written in the form 
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e 1 


r 


i +p{ 



\ s+ p ) . 


( 97 ) 


Evaluating this transfer function with the final value theorem yields a measure of 
steady-state error as the ratio of the magnitude of the steady-state error to the magnitude 
of the reference input signal. 

First, make the substitutions 


a 3 ~ 1 > 


( 98 ) 


a 


2 



( 99 ) 


and 


c k 

a, = — p + — , 
m m 



The transfer function from r to e can be written in the form 


F(s) = 


a 3 s 3 +a 2 s 2 -t-aj5-i-a 0 
ci 3 s^ + d 2 + d 0 + K . 


( 100 ) 


( 101 ) 


( 102 ) 
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To evaluate the steady-state error for a unit step input from r, the term (1/s) must be 
added to account for the step input. Thus, 


e(s ) = F(s)r(s) = 


a 3 s 3 + a 2 s 2 + a x s + a 0 1 
a 3 s + a 2 s 2 +a 1 s + a 0 +K s . 


(103) 



where e ss is the steady-state error with a reference input of a unit step. It is seen that for 
p=0, there is no steady-state error. Also, the steady-state error increases as p increases for 
a given value of K. This is verified through time domain simulation for the case where 
c=10, k=l, m=l, K=l, and p=0.1. Based on Equation ( 104 ), a steady-state error of 
0.091 is expected for a unit step input in r. The step response is shown in Figure 3-23 for 
the system with a pure integrator and the system with a modified integrator. 
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Figure 3-23. Step Response Comparison for Pure and Modified Integrator 


It is seen from Figure 3-23 that the system with the pure integrator does not have 
steady-state error, while the system with the modified integrator has steady-state error of 
0.091 as expected. However the system with the pure integrator sacrifices system 
robustness when compared to this system with the pure integrator. It is demonstrated in 
this section that the system with the modified integrator is significantly more robust to 
error in the sensor measurement. Sensor error is introduced as shown in the block 


diagram in Figure 3-24. 
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Figure 3-24. System with Sensor Error 


Applying a sensor gain of 15, the systems are compared. The response of the 
variable y is shown in Figure 3-25 for each system when r is a unit step input. 


Comparison of Pure/Modified Integral Control 
Sensor Gain = 16 



0 10 20 30 40 50 60 70 80 

Time 


Figure 3-25. Comparison of Pure/Modified Integral Control 
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While neither system follows the step input very well, it is clear that the system 
with pure integral control is unstable, while the system with modified integral control is 
stable. In fact, the gain margin for the system with pure integral control is 10 (20 dB), 
while the gain margin for the system with modified integral control is 20.1 (26 dB). It is 
notable that the system with the modified integrator has a gain margin twice as large as 
the system with pure integral control while sacrificing steady-state error of 9.1%. 

In the case of the centrifuge rotor, the steady-state error analysis is modified to 
treat the MIMO case. If the step input r = r amp * (1/s), then e = F * (1/s) * r amp . Using the 
induced 2-norm (denoted as i2) to define performance, the steady-state error can be 
expressed as 


klL s l F < s =°)|| n H 2 , < 105 > 

where the symbol F denotes the closed-loop transfer function from r (or the disturbance 
in the case of the centrifuge rotor) to e (the translation in the case of the centrifuge rotor). 
This equation is valid for the case where r is a signal with finite area under the curve. In 
the case where r is a step input, the relationship is written using the induced 1-norm 
(denoted as il) [32] as 


^IL^| F ^ = °)|L r an,p 


(106) 


The steady-state error of the centrifuge rotor with a controller is analyzed in Section 


4.3. 1.2. 



66 


3.6 A Constrained Optimization Approach to Control Design 

To design the ABC a procedure is proposed to combine both time domain and 
frequency domain constraints [33] [34], The ABC design problem is first converted into a 
constrained optimization problem where the translation of the rotor due to a fixed worst- 
case imbalance during “spin-up” is minimized subject to problem constraints. A 
constrained optimization approach is then applied to solve the problem. 

A flow chart of the process used in a robust constrained optimization approach to 
control design is shown in Figure 3-26. This process is applied to the mathematical 
formulation of the control problem presented in Section 4.2 to generate an ABC. The 
optimization problem is solved using a gradient-descent method [38] [39] [40]. 



Figure 3-26. Flow of Robust Constrained Optimization Approach to Control Design 
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This approach has also been used to design control systems for the International 
Space Station Momentum Manager, an active vibration isolation system, and the Control 
Momentum Gyroscope Flex Filter for the International Space Station [35] [36] [37]. A 
major advantage of this approach is the capacity to deal with both time domain 
constraints (such as rotor translation constraints) and frequency domain constraints (such 
as robust stability specifications) simultaneously. 
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4 Solution Method 

To solve the ABC design problem, a robust constrained optimization approach to 
control design is implemented. The problem to be solved is re-stated more specifically in 
Section 4.1 and Section 4.2. A three stage design procedure is implemented and the 
intermediate design results are presented in Section 4.3. 

4.1 Problem Setup 

It is assumed that the following specific constraints apply to the ABC: 

I. limit the peak translation (Jg 2 + 77 2 ) of the rotor < 3 . 5 mm for all rodent 
disturbances, 

II. stability is required over the range 0 < co < aw x > 

III. be robust to uncertainty in: 

a. spin rate (to) of 100 %, 

b. rotor mass (M) of 60%, and 

c. rotor inertia (Id,I z ) of 45%, 

IV. limit commanded balancer mass speed < lmm/sec, and 

V. limit steady-state rotor translation < 0.1mm for a fixed disturbance 
(rodents are stationary). 

Note that the tilt of the rotor is not directly constrained since it contributes to the 
translation of the rotor (which is constrained). The above formulation is a specific 
statement of the problem to be solved from the more general problem statement presented 


in Section 2.3. 



4.2 Mathematical Problem Formulation 

The mathematical formulation of the optimization problem can then be posed as: 


Minimize .. ,, 

max i1 2 

(Objective Function) 

(107) 

subject to 

MA-BK)< 0 

(Nominal Stability) , 

(108) 

max||r|| 2 < 3.5 mm 

(Nominal Performance) , 

(109) 

K^)|<0.1mm 

(Nominal Performance) , 

(110) 

ju rs {GU<»)) < 1 

(Robust Stability) , 

(111) 

ju rp {G(joj))<l 

(Robust Performance) , 

(112) 

and 



db < i mm 
dt sec 

(Actuator Limits) , 

(113) 


where r(t) = [£(t) Tj{t)J and b = \{ bl r\ bl £ b3 f b4 ] T . The symbols A and B denote 

the state-space matrices of the coupled plant dynamics of Equations ( 52 ) and ( 53 ). The 
symbol G is defined as the closed-loop centrifuge rotor system as shown in Figure 3-20 
where the uncertainty is in the spin rate, rotor mass, and rotor inertia as defined in Section 
4.1. In this formulation, the peak of the translational response of the rotor is minimized. 
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4.3 Three Stage Solution Procedure 

The controller design task is broken into three separate stages. This is done to 
reduce the number of degrees-of-freedom and constraints in the optimization problem. 
This approach is called “constraint relaxation” in the literature [ 41 ]. Since the equations 
of motion of the' plant can be easily decoupled (see Section 2 . 1 ), it is possible to design a 
controller for the translational equations of motion (Equation ( 54 )) and a controller for 
the tilting equations of motion (Equation ( 55 )) separately, making up the first two stages 
of the design. The third stage then involves designing a controller for the coupled 
dynamics (Equations ( 52 ) and ( 53 )), using the controllers resulting from stage 1 and 
stage 2 as a starting point. Stage 3 is used to fill in the cross terms of the controller as 
shown in Equation ( 117 ). 

The final controller has the structure 

b = Kx ( 114 ) 

where 

b = [Cbl tfbl £b3 *?bA ’ ( 115 ) 

x = [c Tf <p ( C i <t>s jC jV<r J*r > (H 6 ) 

and 

K tran 0 , K crassly 1 

K= 2x6 ~ 2x6 | ( 117 ) 

_K _cross2 2x6 K _tilt 2x6 J 
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K_tran is designed in stage 1 , K_ti.lt is designed in stage 2, and K_crossl and K_cross2 
are designed in stage 3. Overall, the controller design has 48 degrees of freedom (the 
control gains of the K matrix), however, breaking the design task into three stages, it is 
re-formulated as three separate optimization problems having 12, 12, and 24 degrees of 
freedom, respectively. This reduction of the problem size allows the optimization routine 
to converge more quickly and reliably. 

4.3.1 Stage 1 


As mentioned in Section 4.3, stage 1 involves the design of K_tran. After stage 1 
is complete, the controller has the structure 


K = 


K _tran 2x6 

®2x6 



( 118 ) 


The objective of stage 1 is to find K_tran such that the peak translation of the 
rotor is mi nimiz ed, while ensuring that the problem constraints are met. The peak 
translation is defined by the maximum translation of the rotor during time domain 
simulation. The scenario simulated includes a spin-up from 0 to the maximum spin rate. 
The maximum spin rate is reached in 1800 seconds and the rotor continues to spin at this 
rate until the simulation ends at 2500 seconds. The disturbance acting on the plant is the 
disturbance resulting from the rodents being stationary in positions causing the largest 
imbalance in the rotor (see Section 2.1). In addition to minimizing the translation, the 
constraints in Equations ( 108 ) — ( 113 ) must be met as well. The mathematical 
formulation is presented in Section 4.3. 1.1. 
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4.3.1.1 Stage 1 Problem Formulation 


Minimize „ „ 

max rll 2 

(Objective Function) 

(119) 

subject to 

MA-BK)< 0 

(Nominal Stability) , 

(120) 

max||r|| 2 < 3.5mm 

(Nominal Performance) , 

(121) 

|K^/)|< 0.1mm 

(Nominal Performance) , 

(122) 

ju rs {G(j( 0 ))<l 

(Robust Stability) , 

(123) 

ju rp {G(ja>))< 1 

(Robust Performance) , 

(124) 

and 



db . mm 

— < 1 

dt sec 

(Actuator Limits) , 

(125) 


where r(t) — \£(t) Tj(t)J and b = \£ hl Tj h2 f . The symbols A and B denote the state- 

space matrices of the full coupled dynamics in Equations ( 52 ) and ( 53 ). The coupled 
dynamic equations are used to check nominal stability to ensure that the controller will 
not be unstable for the coupled system when stage 3 of the design process is reached. 
The robust stability and robust performance measures must be met for the range 0 < co < 
©max and 60% rotor mass uncertainty from the nominal value. G is defined as the closed- 
loop centrifuge rotor system as shown in Figure 3-20 where the uncertainty is in the spin 


rate and rotor mass as defined above. 
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4.3.1. 2 Stage 1 Design Loop Results 

The constrained optimization problem was solved using the Matlab function 
fmincon [38]. The initial guess to start the optimization process is K_tran = 02x6- The 
optimization took over 1000 iterations to yield the final solution which is presented in 
this section. 

From Figure 4-1, it is apparent that the closed-loop system achieves nominal 
stability since all of the poles lie to the left of the j co-axis. The poles that appear to lie 
near the j co-axis are the result of the modified integrators that were appended to the 
system. Since the modified integrators add poles at s = -0.0003 to the open-loop system, 
these poles appear to be quite near the imaginary-axis, but are in the left half plane 
(LHP). 


Pole-Zero Map 



Figure 4-1. Stage 1 - Nominal Stability 
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Figure 4-2 shows the translation of the rotor when simulated under the conditions 
described in Section 4.3.1. The peak translation of the rotor is well under the maximum 
allowable translation 3.5mm. In fact, the peak translation is 2.52mm. Further, the 
steady-state translation of the rotor is less than 0.10mm satisfying the requirement for 
steady-state error (the constraint imposed on the final value of rotor translation). 



Figure 4-2. Stage 1 - Rotor Translation 
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The steady-state error that is present in the response is due to the modified 
integrators being used in the controller (see Section 3.5). Recall that in the MIMO case 
steady-state error for a step input can be expressed as 

|K,|L<||F( S = 0)y|r„,|^ (126) 


In the case of the centrifuge rotor, the symbol F denotes the closed-loop transfer 
function from the rodent disturbance to rotor translation. The closed-loop centrifuge 
rotor system for stage 1 yields 



-0.00002 

-0.00023 

0.00021 

-0.00003 


0.3 

0.3 


oo 


= 0.000078 m , 


(127) 


where the maximum size of the elements in the vector r is 0.3. Since the maximum size 
of the steady-state error is 0.000078 m in each direction (£ and rj), the upper bound on 
the translation is (0.000078 2 + 0.000078 2 ) 172 = 0.0001 103m. From Figure 4-2, it is seen 
that this relationship provides a good upper bound. The actual steady-state error is 
slightly less than 0.0001m. 
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Figure 4-3 indicates that the closed-loop system maintains stability for 60% mass 
uncertainty over the range 0 < (0 < GW since the upper bound of p < 1 at all frequencies. 
The interval for G) was divided into several intervals to evaluate robust stability as shown 
in Section 3.4.3. The plot of robust stability presented is the interval with the highest 
peak for the upper bound of p. 

Upper Bound on |i - Robust Stability 



Figure 4-3. Stage 1 - Robust Stability 
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The transfer function evaluated to analyze robust stability is the transfer function 
from ua to yA shown in Figure 4-4. The uncertainty A is normalized to 1, and it 
represents the uncertainty present in the rotor mass and rotor spin rate. The symbols D 
and D' 1 are as specified in Section 3.4 and represent scaling matrices used by the 
structured singular value. 



Ya 


Figure 4-4. Block Diagram for Evaluating Robust Stability 
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Figure 4-5 indicates that the closed-loop system achieves the performance 
objective for 60% mass uncertainty over the range 0 < © < gw since the upper bound of 
ft < 1 at all frequencies. The interval for © was again divided into several intervals. The 
plot of robust performance presented is for the © interval with the highest peak for (I. 


Upper Bound on |i -Robust Performance 
60% Mass Uncertainty 



The robust performance objective was to limit rotor translation < 3.5mm for a 
disturbance of maximum amplitude (as in the time domain simulation) and frequency 
content from DC to 0.0001 Hz (0.000628 rad/sec). This performance objective is 
formulated by appending weighting functions to the unweighted plant as in Section 3.4.2. 
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The weighting functions are appended to the plant as shown in Figure 4-6 to 
normalize the performance objective as stated above from d to e to a performance 
objective from d’ to e’ as shown in Figure 4-6. The transfer function evaluated for robust 

performance is [ua d’] T -> [yA e’] T 



Figure 4-6. Weighting Functions for Robust Performance 
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The weighting functions used to normalize the performance objective are 
presented in Figure 4-7 and Figure 4-8 where 


W'. = 

in 


W„ 


i 1 

w. 'I 

rr in J 


and 


W' 


out 



(128) 


(129) 


The weighting function rolls off near 0.000628 rad/sec since the performance 
objective is to reject disturbances with frequency content less than 0.0001 Hz (0.000628 
rad/sec). The magnitude of the weighting function is 0.3 is since the largest expected 
disturbance has a magnitude of 0.3. 


Stage 1 Weighting Functions - W jn 



Figure 4-7. Stage 1 Weighting Functions - Input Weight 
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Figure 4-8 shows the weighting function used to normalize the output to 1. Since 
the largest allowable rotor translation is 3.5mm or 0.0035m at all frequencies, the 
weighting function has the magnitude 285.7. 


Stage 1 Weighting Functions - W out 



Figure 4-8. Stage 1 Weighting Functions - Output Weight 
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Figure 4-9 shows that in the time domain simulation the actuator saturation 
constraints of lmm/s (O.OOlm/s) are not exceeded. By avoiding actuator saturation, the 
linearity of the system is preserved. This is important because linear analysis results 
(such as robust stability and robust performance) are not valid if the actuator constraints 
are saturated since a system with saturation is not linear. 



i i i i i i 

0 500 1000 1500 2000 2500 


Time(s) 


Figure 4-9. Stage 1 - Actuator Constraints 


It is seen that all of the constraints are met in stage 1 and the peak rotor translation 


is minimized. 
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4.3.2 Stage 2 


As stated in Section 4.3, stage 2 involves the design of K_tilt. After stage 2 is 
complete, the controller has the form 


K = 


K _tran 2x6 

®2x6 


^2x6 ,1 

K _ tilt 2x6 J 


(130) 


The objective of stage 2 is to find K_tilt such that the peak tilt of the rotor is 
minimized, while ensuring that the problem constraints are met. The peak tilt is defined 
by the maximum tilt of the rotor during time domain simulation. The scenario simulated 
is the same as that in stage 1. In addition to minimizing the tilt, the stability, robustness, 
steady-state error, and actuator constraints must be met as well. The mathematical 
formulation is presented in Section 4.3.2. 1 . 
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4.3.2.1 Stage 2 Problem Formulation 


Minimize „ „ 

max|4 

(Objective Function) 

(131) 

subject to 

A(A - BK) < 0 

(Nominal Stability) , 

(132) 

max||0|| < 0.6 mrad 

(Nominal Performance) , 

(133) 

\</>(t fma ,)\<0A m rad 

(Nominal Performance) , 

(134) 

M„{CHja»)< 1 

(Robust Stability) , 

(135) 

ju rp {G(jO)))< 1 

(Robust Performance) , 

(136) 

and 



db ,mm 
— <1 — 
dt sec 

(Actuator Limits) , 

(137) 


where </>(t) = [fy(t) and b =[£ M ] r . The symbols A and B denote the state- 

space matrices of the full coupled dynamics in Equations ( 52 ) and ( 53 ). The robust 
stability and robust performance measures must be met for the range of 0 < co < ow and 
30% independent inertia uncertainty in Id and I z . This inertia uncertainty will be 
increased to 45% in stage 3. G is defined as the closed-loop centrifuge rotor system as 
shown in Figure 3-20 where the uncertainty is in the spin rate and rotor inertia as defined 
above. The tilt of the rotor is not constrained in Equations ( 108 ) — ( 113 ) since in the 
coupled system it contributes a small amount to the translation. However, that is not the 
case in the uncoupled system used for stage 2 of the controller design, so the tilt is 
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constrained to be < 0.6 mrad (which would contribute 0.35mm of translation in the 
coupled system or 10% of the maximum allowable translation). 

4.3.2.2 Stage 2 Design Loop Results 

The constrained optimization problem was solved using the Matlab function 
fmincon [38]. The initial guess to start the optimization process is K_tilt = 02x6- The 
optimization took over 1000 iterations to yield the final solution which will be presented 
in this section. 

From Figure 4-10, it is apparent that the closed-loop system achieves nominal 
stability since all of the poles lie to the left of the joo-axis. As in Figure 4-1 the poles that 
appear to lie near the j co-axis are in the LHP and are the result of the modified integrators 
that were appended to the system. 
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Figure 4-10. Stage 2 - Nominal Stability 
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Figure 4-11 shows the tilt of the rotor when simulated under the conditions 
described in Section 4.3.1. The peak tilt of the rotor is less than the maximum allowable 
tilt of 0.6mrad. In fact, the peak tilt is 0.54mrad. Further, the steady-state tilt of the rotor 
is less than O.lOmrad, satisfying the requirement for steady-state error (the constraint 
imposed on the final value of rotor tilt). 



Figure 4-11. Stage 2 - Rotor Tilt 
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Figure 4-12 indicates that the closed-loop system maintains stability for 30% 
inertia uncertainty over the range 0 < co < ©max since the upper bound of p < 1 at all 
frequencies. Again, the interval for © was divided into several intervals to evaluate 
robust stability as shown in Section 3.4.3. The plot of robust stability presented is the 
interval with the highest peak for the upper bound of p. The transfer function evaluated 
to analyze robust stability is the same as that in Figure 4-4, but the uncertainty block 
A represents uncertainty in the rotor inertia and rotor spin rate. The plant is the tilting 
dynamics in Equation ( 55 ) and the controller is K_tilt. 


Upper Bound on |i - Robust Stability 
30% Inertia Uncertainty 



Figure 4-12. Stage 2 - Robust Stability 
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Figure 4-13 shows that the closed-loop system achieves the performance objective 
for 30% inertia uncertainty over the range 0 < (0 < ow since the upper bound of |l < 1 at 
all frequencies. The interval for to was again divided into several intervals. The plot of 
robust performance presented is for the to interval with the highest peak for the upper 
bound of \x. 


Upper Bound on p - Robust Performance 
30% Inertia Uncertainty 



The robust performance objective in stage 2 was to limit the tilt of the rotor < 
0.6mrad for a disturbance of maximum amplitude (as in the time domain simulation) and 
frequency content from DC to 0.0001 Hz (0.000628 rad/sec). This performance objective 
is formulated by appending weighting functions to the unweighted plant. 
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The weighting functions are appended to the plant as shown in Figure 4-6. Again 
the transfer function evaluated is [ua d’] T -* [yA e’] T from Figure 4-6. The weighting 

functions used to normalize the performance objective for stage 2 are presented in Figure 
4-14 and Figure 4-15 where 


and 


W'. = 

r in 


W in 1 

in I i 

W. I 

in J 


W' = 

rr out 




W M 

out J 


(138) 


(139) 


The weighting function for the input to the system rolls off near 0.000628 rad/sec 
since the performance objective is to reject disturbances with frequency content less than 
that. The magnitude of the weighting functions at low frequencies is 0.00033 since that 
is the magnitude of the largest expected disturbance in the frequency range. 



Figure 4-14. Stage 2 Weighting Functions - Input Weight 
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Figure 4-15 shows the weighting function used to normalize the output to 1. 
Since the largest allowable rotor tilt is 0.6mrad or 0.0006rad at all frequencies, the 
weighting function has the magnitude 1666.7. 

Stage 2 Weighting Functions - W 



Frequency (rad/sec) 

Figure 4-15. Stage 2 Weighting Function - Output Weight 
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Figure 4-16 shows that in the time domain simulation, the actuator saturation 
constraints of 1 mm/sec or 0.001 m/sec are not exceeded. In fact, the mass rates are far 
from the saturation limits. The lower plot is a closer view of the mass rates. By avoiding 
actuator saturation, the linearity of the system is preserved. 
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Figure 4-16. Stage 2 - Actuator Constraints 


It is seen that all of the constraints are met in stage 2 and the peak rotor tilt is 


minimized. 
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4.3.3 Stage 3 

As mentioned in Section 4.3, stage 3 involves the design of K_crossl and 
K_cross2. These sub-matrices are designed simultaneously allowing for 24 degrees of 
freedom (12 for each K_cross matrix) in the stage 3 problem formulation After stage 3 
is complete, the controller has the form 


K _tran 2x6 K _cross\ 2x6 ~\ 
K _cross2 2x6 K _ tilt 2x6 J 


(140) 


The objective of stage 3 is to find K_crossl and K_cross2 such that the peak 
translation of the rotor is minimized, while ensuring that the problem constraints are met. 
Further, it is desirable to increase the amount of inertia uncertainty that the system can 
tolerate. The requirement for inertia uncertainty is increased to 45% uncertainty in both 
Id and I z . In addition to minimizing the translation, the constraints in Equations ( 108 ) - 
( 113 ) must be met as well. The mathematical formulation is presented in Section 
4.3.3. 1. This formulation is the same as formulation in Section 4.2 since stage 3 is the 
last stage in the design process and all problem constraints must be met for the system 
with coupled dynamics. 
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4.3.3.1 Stage 3 Problem Formulation 


Minimize „ „ 

max|r| 2 

(Objective Function) 

(141) 

subject to 
X(A-BK)<0 

(Nominal Stability) , 

(142) 

max||r|| 2 < 3.5 ram 

(Nominal Performance) , 

(143) 

|K^„ a/ )|<0.1mm 

(Nominal Performance) , 

(144) 

ju rs {G (ja»)<l 

(Robust Stability) , 

(145) 

Prp (G(jCt))) < 1 

(Robust Performance) , 

(146) 

and 



db . mm 

— <1 

dt sec 

(Actuator Limits) , 

(147) 


where r(t) = [£(t) ^(ofand b = [C bl J] b2 g b3 % M J . The symbols A and B denote 
the state-space matrices of the full coupled dynamics in Equations ( 52 ) and ( 53 ). The 
robust stability and robust performance measures must be met for the range of 0 < (0 < 
tOmax, 45% of rotor inertia uncertainty in both I z and Id, and 60% rotor mass uncertainty 
from the nominal value. G is defined as the closed-loop centrifuge rotor system as shown 
in Figure 3-20 where the uncertainty is in the spin rate, rotor mass, and rotor inertia as 


defined above. 
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4.3.3.2 Stage 3 Design Loop Results 

The constrained optimization problem was solved using the Matlab function 
fmincon [38]. The initial guess to start the optimization process is K_crossl = K_cross2 
= 02x6- The optimization took over 1000 iterations to yield the final solution which will 
be presented in this section. 

From Figure 4-17, it is apparent that the closed-loop system achieves nominal 
stability since all of the poles lie to the left of the jco-axis. As in Figure 4-1 and Figure 
4-10 the poles that appear to lie near the jco-axis are in the LHP and are the result of the 
modified integrators that were appended to the system. 


Pole-Zero Map 



Figure 4-17. Stage 3 - Nominal Stability 




95 


Figure 4-18 shows the translation of the rotor when simulated under the 
conditions described in Section 4.3.1. The peak translation of the rotor is well under the 
maximum allowable translation of 3.5mm. In fact, the peak translation is 2.54mm. It is 
seen that the peak translation is slightly higher than the final result of stage 1 . This is due 
to the fact that the coupled dynamic equations of motion include the contribution from 
the tilt in the rotor translation. Further, the steady-state translation of the rotor is less than 
0.10mm satisfying the requirement for steady-state error (the constraint imposed on the 
final value of rotor translation). 



Figure 4-18. Stage 3 - Rotor Translation 
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Figure 4-19 indicates that the closed-loop system maintains stability for 60% rotor 
mass uncertainty and 45% rotor inertia uncertainty over the range 0 < (0 < ow since the 
upper bound of |i < 1 at all frequencies. The interval for to was again divided into several 
intervals to evaluate robust. The plot of robust stability presented is the interval with the 
highest peak for the upper bound of JJ,. The transfer function evaluated is the same as that 
evaluated in stage 1 and stage 2 of the design process, but the plant includes the full 
coupled dynamics and the uncertainty is in the rotor mass, rotor inertia, and rotor spin 
rate. 


Upper Bound on |i - Robust Stability 
60% Mass/45% Inertia Uncertainty 



Figure 4-19. Stage 3 - Robust Stability 
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Figure 4-20 shows that the closed-loop system achieves the performance objective 
for 60% rotor mass uncertainty and 45% rotor inertia over the range 0 < to < ow since 
the upper bound of p < 1 at all frequencies. Again, the interval for co was again divided 
into several intervals. The plot of robust performance presented is for the co interval with 
the highest peak for the upper bound on p. The transfer function evaluated is the same as 
that in stage 1 and stage 2 of the design process, but the plant includes the full coupled 
dynamics and the uncertainty is in the rotor mass, rotor inertia, and rotor spin rate. 


Upper Bound on |i - Robust Performance 
60% Mass/45% Inertia Uncertainty 



Figure 4-20. Stage 3 - Robust Performance 
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The robust performance objective was to limit the translation of the rotor < 
3.5mm for a disturbance of maximum amplitude (as in the time domain simulation) and 
frequency content from DC to 0.0001 Hz (0.000628 rad/sec). This performance objective 
is formulated by appending weighting functions to the unweighted plant. 

The weighting functions are appended to the plant as shown in Figure 4-6 in order 
to normalize the performance objective as stated above from d to e to a performance 
objective from d’ to e’. 

The weighting functions used to normalize the performance objective for stage 3 
are a combination of the weighting functions for stage 1 and stage 2. Specifically, 


W, 


in stage! 


w, 


in stage 1 


w ; 


in stage 1 


W, 


in stage 2 


W, 


in stage 2 


(148) 


and 


W. 


out stage 3 




out stage 1 


w. 


out stage 1 


W, 


out stage 2 


w. 


out stage 2 


(149) 
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Figure 4-21 shows that in the time domain simulation the actuator saturation 
constraints of lmm/sec or O.OOlm/sec are not exceeded. By avoiding actuator saturation 
the linearity of the system is preserved. 


x 1 g 3 Commanded Mass Rate - Constraints 



Figure 4-21. Stage 3 - Actuator Constraints 


It is seen that all of the problem constraints specified in Equations ( 108 )-( 113) 
are met in stage 3 and the peak translational response of the rotor is minimized resulting 
in the final design of the ABC. 
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5 Results 

In order to further examine of the effectiveness of the ABC design presented in 
Chapter 4, analysis was performed using frequency domain techniques to further analyze 
stability and time domain simulation (including Monte Carlo simulation) to further 
analyze performance. Much of the needed analysis was done during the design phase in 
Chapter 4 where nominal stability, nominal performance, robust stability, robust 
performance, and actuator saturation are analyzed for the linear system in the design 
phase. In this chapter, both stability margin analysis and nonlinear time domain analysis 
are performed to further verify the ABC design. 

5.1 Stability Margin 

Frequency domain stability and performance results are seen in Figure 4-19 and 
Figure 4-20 in Chapter 4 for the final ABC. These results guarantee the stability and 
performance of the closed-loop system for the range of possible spin rates and for 
parametric variations of 60% in the mass of the rotor (M) and 45% in the transverse 
moment of inertia (Id) and polar moment of inertia (I z ) values. These requirements were 
the primary stability constraints of the control design. This section provides further 
verification of the stability of the system. 

As an additional measure of the stability of the closed-loop system, gain margin 
and phase margin are calculated. Gain margin refers to the amount of gain that could be 
multiplied by the signal at the point where the loop break is shown while maintaining 
closed-loop stability. Phase margin refers to the amount phase that could be added to the 
signal at the point where the loop break is shown while maintaining closed-loop stability. 
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The control system is analyzed one loop at a time. The loops are shown in Figure 
5-1 and Figure 5-2. First, the stability margins are evaluated at the input to the plant (in 
the actuator signal). Next, the stability margins are evaluated at the output of the plant (in 
the sensor signal). 

The stability margins at the input of the plant as shown in Figure 5-1 are 
presented in Table 5-1. All stability margins are computed for the plant with a nominal 
spin rate, C0ma X /2. 


Loop Break 



Figure 5-1. Loop Break for SISO Stability Margins at Plant Input 


Traditional guidelines for good stability are gain/phase margins of 6dB/30deg or 
better. From the results in Table 5-1, it is apparent that the gain margins far exceed 6dB, 
with the smallest gain margin in the ^3 -*■ ^>3 loop being 18.5dB. The smallest phase 

margin, 28.9 deg occurs in the r|b 2 - T|b 2 loop. While this is slightly less than 30 deg 


phase margin, it is acceptable. 
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Gain Margin(dB) 

Phase Margin(deg) 

Chi Chi 

34.7 

37.5 

f]b2 -*■ ^2 

31.0 

28.9 

Cb3 -*■ Cb3 

18.5 

102.7 

Cb4 -*■ Ch4 

Inf 

107.8 


Table 5-1. SISO Stability Margins at Plant Input for Nominal Spin Rate 


The stability margins at the output of the plant as shown in Figure 5-2 are 
presented in Table 5-2. All stability margins are computed for the plant with a nominal 
spin rate, GW^. 



Figure 5-2. Loop Break for SISO Stability Margins at Plant Output 




103 


The results in Table 5-2 indicate that the gain margins at the output of the plant 
far exceed the traditional guideline of at least 6dB. Similar to the phase margin in the 
T lb 2 -*- "Hba loop at the input of the plant, the T|_int - T)_int loop at the output of the plant has 

a phase margin of approximately 28 deg. Again, while this is slightly less than 30 deg 
phase margin, it is acceptable. The notations _dot and _int denote the derivative and 
integral of the associated variable. 



Gain Margin(dB) 

Phase Margin(deg) 

c 

55.2 

Inf 

T1 — T| 

44.9 

Inf 

<te 

30.5 

Inf 

-*• <(>n 

45.3 

Inf 

£_dot -» £_dot 

43.4 

Inf 

r)_dot -*■ r]_dot 

37.9 

Inf 

<()f;_dot -» (j)^_dot 

Inf 

Inf 

(|) ri _dot -> (j>n_dot 

21.1 

Inf 

^_int -*■ ^_int 

36.0 

37.5 

r|_int -*• T|_int 

35.3 

27.9 

(J>^_int -»• (he. int 

37.2 

83.1 

<h n _int -> (|)r|_int 

43.6 

100.2 


Table 5-2. SISO Stability Margins at Plant Output for Nominal Spin Rate 


Although gain/phase margins were not constrained during the design process, the 
resulting system has good gain/phase margins. This is primarily a result of two 
constraints that were in the design process. First, the design process required that the 
system be robustly stable to large parametric uncertainty in the rotor mass and the rotor 
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inertia. Since changing these parameters changes the gain of the plant, the controller had 
to be robust to gain uncertainty, yielding acceptable gain margin. Second, the constraints 
on the speed of the balancer masses necessitated that the controller gains stayed relatively 
small. The gains in the controller gain matrix could not be large values since such values 
would cause the controller to command the masses to move at a speed greater that the 
saturation speed. Having a low-gain controller corresponds to good stability margins 
much like a system with a zero gain controller (no control) would have infinite stability 
margins. However, not having a controller is not a reasonable solution since the 
controller is needed to meet time domain performance requirements (see Figure 5-4). 

5.2 Time Domain Nonlinear Simulation 

Time domain performance is analyzed using a simulation as shown in Figure 5-3. 
The simulation includes nonlinearities in the coordinate transformation, actuator 
constraints (mass rate saturation), and time delays. Also, the plant dynamic equations of 
motion are in the fixed (x,y) frame and the sensor measurements are rotated into the 
rotating frame via the coordinate transformation following the plant model. 
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Figure 5-3. Time Domain Simulation 

The scenario simulated in this section includes a spin-up of the rotor from a spin 
rate of 0 to the maximum spin rate. The maximum spin rate is reached in 1800 seconds 
and the rotor continues to spin at this rate until the simulation ends at 2500 seconds. The 
disturbance acting on the plant is the disturbance resulting from the rodents being 
stationary in positions causing the largest imbalance in the rotor as described in Section 
2 . 1 . 

First, the performance of the ABC is analyzed by comparing the rotor translation 
for the case with no controller and the case where the controller designed in Chapter 4 is 
implemented. The results of these simulations are shown in Figure 5-4. 




106 



Time(s) 


Figure 5-4. Time Domain Comparison to No Control 


It is clear that control is needed to satisfy the requirement of limiting peak rotor 
translation to be less than 3.5mm. In fact, with the ABC the peak rotor translation is 
reduced from approximately 7.0mm to approximately 2.5mm. 

Further, the effect of the specified variations in the rotor mass and inertia was 
analyzed through time domain simulation. The scenario simulated was the same as that 
described above that was used to compare the rotor translation with the ABC to a case 
with no control. In order to analyze the effect of parametric variation on the performance 
of the system, the comer cases of the allowed parametric variations were tested. For 
instance, the mass of the rotor was set at the nominal value, the nominal value + 60%, 
and the nominal value - 60%. Similarly, the inertias of the rotor (Id,I z ) are allowed to 
vary independently between the nominal value, the nominal value + 45%, and the 
nominal value - 45%. Using these comer cases, 27 different parameter combinations 
were tested in simulation. The combinations are listed in Table 5-3. Figure 5-5 shows 
the response of all 27 cases tested. The peak rotor translation of the 27 cases varied 
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between 2.07mm and 2.87mm, all below the 3.5mm requirement. Also, the translation of 
the rotor at steady-state varied between 0.06mm and 0.10mm. This meets the 
requirement for steady-state error of 0. 10mm or less. 


Rotor Translation 


x in ' 3 Parametric Variations 



From Table 5-3 it is seen that the combination of parameter values corresponding 
the largest peak rotor translation and the largest steady-state error are when the mass of 
the rotor is +60% from nominal, the transverse moment of inertia of the rotor is +45% 
from nominal, and the polar moment of inertia of the rotor is -45% from nominal. This 
case is highlighted in Table 5-3. 
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M (kg) 

Id (kg- 

Iz (kg- 

Peak Translation 

Steady-State 

m A 2) 

m A 2) 

(mm) 

Translation (mm) 

nom 

nom 

nom 

2.53 

0.093 

-60% 

nom 

nom 

2.20 

0.093 

+60% 

nom 

nom 

2.70 

0.076 

nom 

-45% 

nom 

2.47 

0.081 

-60% 

-45% 

nom 

2.15 

0.087 

+60% 

-45% 

nom 

2.64 

0.062 

nom 

+45% 

nom 

2.59 

0.098 

-60% 

+45% 

nom 

2.25 

0.098 

+60% 

+45% 

nom 

2.76 

0.084 

nom 

nom 

-45% 

2.63 

0.098 

-60% 

nom 

-45% 

2.28 

0.098 

+60% 

nom 

-45% 

2.80 

0.087 

nom 

-45% 

-45% 

2.57 

0.098 

-60% 

-45% 

-45% 

2.22 

0.096 

+60% 

-45% 

-45% 

2.73 

0.081 

nom 

+45% 

-45% 

2.69 

0.099 

-60% 

+45% 

-45% 

2.33 

0.099 

+60% 

+45% 

-45% 

2.87 


nom 

nom 

+45% 

2.44 

0.081 

-60% 

nom 

+45% 

2.12 

0.084 

+60% 

nom 

+45% 

2.65 

0.060 

nom 

-45% 

+45% 

2.39 

0.090 

-60% 

-45% 

+45% 

2.07 

0.086 

+60% 

-45% 

+45% 

2.54 

0.085 

nom 

+45% 

+45% 

2.50 

0.087 

-60% 

+45% 

+45% 

2.17 

0.091 

+60% 

+45% 

+45% 

2.66 

0.066 


Table 5-3. Results of Parametric Variations 


The preceding time domain analysis verifies two things. First, it is apparent that 
without control the system will not meet the performance requirements of limiting rotor 
translation to a peak of less than 3.5mm or a steady-state value of 0.10mm. This verifies 
the need for an ABC to improve performance. Second, parametric variations in time 
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domain simulation verify that the performance of the system is acceptable subject to the 
allowable parametric uncertainty. The system performance is robust to uncertainty in the 
plant. The robustness is further verified through Monte Carlo Simulation. 

5.3 Monte Carlo Simulation 

The performance of the ABC system is further analyzed through Monte Carlo 
simulation. The Monte Carlo method, a numerical method of solving problems through 
the simulation of random variables and processes [42],[43], has been used to gain insight 
into many complex problems in engineering and other fields. The method allows for the 
repeated simulation of a given system with an unknown or uncertain input variable. 
Repeated random selection of this unknown variable allows one to study how a particular 
system’s outputs are distributed given a random input. 

Five separate Monte Carlo studies were performed to analyze the performance 
and robustness of the ABC system subject to uncertainty in the following variables: 
rodent motion within the habitats, mass and inertia of the rotor, actuator uncertainty, and 
sensor uncertainty. 

The Monte Carlo studies were performed under the conditions shown in Table 
5-4. 0. 1 m/s is assumed to be the maximum velocity of a single rodent. 
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Test 
Case • 

Random 

Input 

Variables 

Distribution 

Disturbance 

Output 

Variable 

Test 
Case 1 

Rodent 

velocity 

Normal, mean = 0 m/s, 
3a = 0. 1 m/s 

Random 
variable due to 
rodent motion 

Peak rotor 
translation 

Rodent 

direction 

Uniform, random unit vector 

Test 
Case 2 

Rotor 

mass 

Normal, mean = nominal, 
3a = 60% variation 

Rodents in 
fixed positions 
for max 
imbalance 

Peak rotor 
translation 

Rotor 

inertia 

Normal, mean = nominal, 
3a = 45% variation 

Test 
Case 3 

Actuator 

gain 

Normal, mean = 1 , 
3a = 75% error 

Rodents in 
fixed positions 
for max 
imbalance 

Peak rotor 
translation 

Test 
Case 4 

Sensor 

gain 

Normal, mean = 1, 
3a = 75% error 

Rodents in 
fixed positions 
for max 
imbalance 

Peak rotor 
translation 

Test 
Case 5 

Rodent 

velocity 

Normal, mean = 0 m/s, 
3a = 0.1 m/s 

Random 
variable due to 
rodent motion 

Peak rotor 
translation 

Rodent 

direction 

Uniform, random unit vector 

Rotor 

mass 

Normal, mean = nominal, 
3a = 60% variation 

Rotor 

inertia 

Normal, mean = nominal, 
3a = 45% variation 

Actuator 

gain 

Normal, mean = 1 , 
3a = 75% error 

Sensor 

gain 

Normal, mean = 1, 
3a = 75% error 


Table 5-4. Monte Carlo Test Cases 


The results of the Monte Carlo study are shown in Figure 5-6 - Figure 5-10. 


Ill 


Test Case 1 - Sample Distribution of Rodent Velocity 



Rodent Velocity (m/s) 


Mean Translation = 1.3mm 
Standard Deviation = 0.1 1 mm 


Test Case 1 - Maximum Rotor Translation 


Performance 

Requirement 


1.5 2 2.5 

Rotor Translation (m) 


Figure 5-6. Test Case 1 Results - Random Rodent Motion 

It is seen in test case 1, that the system meets the performance requirement for all 
1000 runs. This result indicates that the ABC is effectively canceling the imbalance 


introduced by rodent motion. 
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Test Case 2 - Sample Distribution of Rotor Mass 

300 1 1 . . . . 1 . — 



400 BOO 1200 1600 2000 2400 2800 

Rotor Mass (kg) 


Fixed 
Disturbance 



Controller 


♦ 


O 


A similar distribution exists 
for the inertias of the Rotor 


Mean Translation = 2.5mm 
Standard Deviation = 0.09mm 


300 r 


Test Case 2 - Maximum Rotor Translation 


performance 

[Requirement 


1.5 2 2.5 

Rotor Translation (m) 


Figure 5-7. Test Case 2 Results - Random Rotor Mass/Inertia 


Similarly, in test case 2, the system meets the performance requirement for all 
1000 runs. This result indicates that the system with the ABC controller is robust to 


variations in the mass and inertia of the rotor. 
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Fixed 

Disturbance 




Plant 



Test Case 3 - Maximum Rotor Translation 



I 0.6% of test cases outside 

Mean Translation = 2.6mm of desired performance 

Standard Deviation = 0.28mm 


Figure 5-8. Test Case 3 Results - Random Actuator Gain 


In test case 3, it is seen that the system meets the performance requirement for 
99.4% of the simulations. It is not surprising that it does not meet the performance 
requirement in every case since a 75% uncertainty in the actuator gain is quite large. 
This allows the balancer mass to be up to 75% away from the position commanded by the 
ABC. A success rate of 99.4% is considered acceptable. 
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Tost Case 4 - Sample Distribution of Sensor Gain 
300 1 . i 



Sensor Gain 


Controller 



Mean Translation = 2.6mm 
Standard Deviation = 0.30mm 


I 
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Test Case 4 - Maximum Rotor Translation 


300 

250 

200 

150 

100 

50 

0 


Performance 



Translation (m) x 10‘ 3 


1 .0% of test cases outside 
of desired performance 


Figure 5-9. Test Case 4 Results - Random Sensor Gain 


In test case 4 it is seen that the system meets the performance requirement for 
99.0% of the simulations. Similar to test case 3, it is not surprising that it does not meet 
the performance requirement in every case since a 75% uncertainty in the sensor gain is 
quite large. This allows the sensor to read a displacement with up to 75% error from the 
actual displacement. A success rate of 99.0% is considered acceptable. 
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A similar distribution exists 
for the inertias of the Rotor 


T*xt C»x* 5 - Sampl* Distribution of Rotor Maxi 


Test Case 5 - Maximum Rotor Translation 


Performance 

Requirement 



0 15 -0 1 -0 05 0 0.05 0 1 0 15 

Rodent Velocity (m/s) 

Test Case 5 - Sample Distribution of Actuator Gain 


Test Case 5 - Sample Distribution of Sensor Gain 



Mean Translation = 1.3mm 
Standard Deviation = 0.15mm 


Controller 


Figure 5-10. Test Case 5 Results - All Random Distributions 


It is seen in test case 5 that the system meets the performance requirement for all 
1000 runs. This result indicates that the ABC is effectively canceling the imbalance 
introduced by rodent motion and is robust to simultaneous perturbations in parameter 
values, actuator gains, and sensor gains. It is not surprising that the mean translation in 
test case 5 is significantly less than that in test cases 2-4. This is due to the fact that the 
disturbance in test cases 2-4 was a fixed disturbance where the rodents are stationary in 
the worst-case positions. The disturbance in test case 5 is a random disturbance 
generated from the random direction and velocity of the rodents. In general, this random 
disturbance will be less severe than the fixed disturbance with maximum imbalance. 





116 


6 Concluding Remarks 

An ABC meeting desired specifications has been developed for the centrifuge 
rotor using a robust constrained optimization approach. This approach is selected 
because constraints in the time domain and frequency domain can easily be handled 
simultaneously with this method. Other ABC design methods do not allow this freedom. 

The effectiveness of the controller has been demonstrated by verifying that the 
constraints are met at the final iteration of the design, checking the stability margins of 
the resulting system, performing selected time domain tests, and testing the robustness of 
the design through Monte Carlo methods. In each case, the performance of the ABC was 
acceptable. 

In Chapter 2 the problem background was introduced. A description of the 
centrifuge rotor system and the ABC system was presented. The equations of motion for 
the centrifuge rotor were derived in both the fixed frame coordinate system and the 
rotating frame coordinate system. The system block diagram and controller structure for 
the ABC were presented. 

In Chapter 3 the technical background for the solution of the problem is 
discussed. The basics of MIMO control and uncertain systems were presented. These 
concepts were used together to analyze system robustness with the structured singular 
value. Steady-state error in control systems was also discussed. Finally the robust 
constrained optimization approach to control design is presented. 

Chapter 4 presented the method used to solve the control problem. This chapter 
described the implementation of the robust constrained optimization approach to control 
design and how it was applied to design an ABC system. The method was applied in a 
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three stage design procedure in order to convert the relatively large optimization problem 
into three smaller sub-problems. The results of the controller design are presented at the 
end of each individual stage and also at the end of the final design. The application of the 
solution method produced a controller meeting the constraints of the problem while 
minimizing the translation of the rotor. 

Chapter 5 presented further analysis of the control design. The system was 
analyzed in both the time domain and frequency domain. Stability margins were 
analyzed in the frequency domain and performance was analyzed in the time domain. 
Frequency domain analysis verified that the stability margins of the system were 
satisfactory. The time domain simulations verified that the controller performance was 
robust to parametric uncertainty. The results provided further indication that the control 
design was acceptable. Also, the results of the Monte Carlo simulations were presented. 
The Monte Carlo results also verified that the control design provided satisfactory 
performance. 

Future work in the design of an ABC for the centrifuge rotor system could include 
the minimization of a different performance measure, such as robust performance or 
robust stability. Such a controller would be able to tolerate more uncertainty, but the 
performance at the nominal point would suffer as a result. Also, a method could be used 
to treat robust performance in the time domain. This would make it possible to consider 
transient effects in the robust performance analysis. 
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